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Abstract 

Recent works have highlighted a double role for the Transforming Growth Factor /? (TGF-p): it inhibits cancer in healthy 
cells and potentiates tumor progression during late stage of tumorigenicity, respectively; therefore it has been termed the 
"Jekyll and Hyde" of cancer or, alternatively, an "excellent servant but a bad master". It remains unclear how this molecule 
could have the two opposite behaviours. In this work, we propose a TGF-p multi scale mathematical model at molecular, 
cellular and tissue scales. The multi scalar behaviours of the TGF-p are described by three coupled models built up together 
which can approximatively be related to distinct microscopic, mesoscopic, and macroscopic scales, respectively. We first 
model the dynamics of TGF-p at the single-cell level by taking into account the intracellular and extracellular balance and 
the autocrine and paracrine behaviour of TGF-p. Then we use the average estimates of the TGF-p from the first model to 
understand its dynamics in a model of duct breast tissue. Although the cellular model and the tissue model describe 
phenomena at different time scales, their cumulative dynamics explain the changes in the role of TGF-p in the progression 
from healthy to pre-tumoral to cancer. We estimate various parameters by using available gene expression datasets. Despite 
the fact that our model does not describe an explicit tissue geometry, it provides quantitative inference on the stage and 
progression of breast cancer tissue invasion that could be compared with epidemiological data in literature. Finally in the 
last model, we investigated the invasion of breast cancer cells in the bone niches and the subsequent disregulation of bone 
remodeling processes. The bone model provides an effective description of the bone dynamics in healthy and early stages 
cancer conditions and offers an evolutionary ecological perspective of the dynamics of the competition between cancer and 
healthy cells. 
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Introduction 

A full systemic understanding of cancer process will benefit from 
investigating cell-tissue interaction. We can observe what happens 
at more or less all scales, from the disease at the whole organism 
down to the molecular level of cancer, and we have good amount 
of experimental data on all levels of biological organization. 
However, putting things together in order to obtain real 
understanding is much more difficult and much less developed. 
A way to build up multi scale models is by using proteins that are: 
1) mutational drivers, meaning the mutation of one of the related 
genes causes the change of the phenotype, 2) able to interact with 
proteins which have intracellular and extracellular effects; hence, 
involving multi-cellular phenomena. Here, we start with the 
consideration that tissue modeling is the missing link between basic 
research and cUnical practice, and we aim at using a modeling 
approach to bridge the cell to tissue scale in health and disease 
(cancer) dynamics. A key player of this multi scale process is TGF- 
P family of cytokines that control numerous cellular responses, 
including proliferation, differentiation, apoptosis and migration. 
TGF-P is always produced as an inactive cytokine that cannot 
bind to its receptor and function unless the latent complex is 
somehow activated. This regulation provides a complex control of 
TGF-p function, thereby ensuring that its potent effects are 



produced in appropriate locations and times. TGF-P interacts 
with cytoskeleton, epitheUal cadherin (E-cad) and integrins 
producing a multi scale mechanobiological effects on tissue [1]. 
Cancer is a multi scale, multifactorial and multi step process [2,3]. 
The cancer cells undergo a cascade of mutations, some of them 
changing the phenotype, to obtain the ability to metastasise, and 
are constantly exposed to signals that induce apoptosis. Acquisition 
of antiapoptotic properties by cancer cells is important for 
metastasis, and recent studies suggest that TGF-P promotes the 
survival of certain types of cancer cells [4,5]. TGF-P both inhibits 
and facilitates tumor progression during early and late stage of 
tumorigenicity, respectively. However, it stiU remains veiled how 
TGF-P plays both contrasting roles [6-8]. Therapies based on 
TGF-P seem promising [9]. Tumorigenesis is in many respects a 
process of disregulated cellular evolution that drives malignant 
cells to acquire several phenotypic hallmarks of cancer, including 
the ability of growing autonomously, disregarding cytostatic 
signals, ignoring apoptotic signals, stimulating angiogenesis, 
invading, metastasising and becoming immortal. In the next 
section, we introduce the role of TGF-P in breast cancer. 

The Ductal Lobular Unit and Breast Cancer 

The terminal ductal lobular unit is the basic functional and 
histopathological unit of the breast, and it has been identified as 



PLOS ONE I www.plosone.org 



1 



February 2014 | Volume 9 | Issue 2 | e88533 



Modeling TGF-p in Cancer Tissue Dynamics 



the site of origin of the most common breast malignancy. The 
ductal carcinoma corresponds to a specific stage of cancer 
development of the mammary parenchyma, Figure (1). Recent 
works showed that TGF-f! is abundantly expressed by highly 
metastatic breast cancer cells and promotes their survival. In 
particular, TGF-p autocrine signaling, in certain breast cancers, 
promotes cell survival via inhibition of apoptotic signaling [10]. 
Major determinants of the "tissue identity" are the cadherins and 
integrins which are adhesion molecules regulating cell-cell and 
cell-matrix interactions. Cells containing a specific cadherin 
subtype tend to cluster together to the exclusion of other types, 
both in cell culture and during development. In vitro and in vivo 
studies have demonstrated the existence of crosstalk between 
integrins and cadherins in cell adhesion and motility. 

Integrins play a key role in activating the latent complex of 
TGF-fi. The binding of integrins to the latent complex induces a 
conformational change that produces active TGF-ji that binds to 
its receptor. Evidence now points to a crucial role for cell 
contraction in the activation of TGF-fi via few types of integrins. 
An intact actin cytoskeleton is required for TGF-fi activation by 
these integrins; activation is greatly reduced by treatment of cells 
with the actin polymerization inhibitor cytochalasin. Additionally, 
treatment of cells with inhibitors of cell contraction, greatiy 
reduces integrin-mediated TGF-j) activation, whereas agents that 
stimulate cell contraction, such as thrombin, angiotensin-II and 
endothelin-1 enhance TGF-j) activation by integrins [1]. 

In breast cancer, the expression of E-cad is a hallmark of a well 
differentiated epithelium that functions to maintain cell-cell 
junctions, thereby inhibiting aberrant cell proliferation and 
migration. The loss of E-cad function via genetic inactivation by 
TGF-fi or via epigenetic silencing is a common characteristic of 
systemicaUy invasive carcinomas. Down-regulated E-cad expres- 
sion is required for the outgrowth of breast cancer cells. Breast 
cancer cells show a typical pattern of dissemination by 1) down- 
regulating E-cad expression or activity; 2) separating cell-cell 
junctions; 3) invading the surrounding tissues; and 4) intravasating 
the vasculature or lymphatic system [1 1]. Recent works show that 
down-regulated E-cad expression induced by TGF-fi was suffi- 
cient to prevent mammary epithelial cell differentiation and 
inst(;ad, produced dense and more spherical cultures that 
underwent metastatic outgrowth [12] [13]. 

Here, we model the molecular and cellular mechanisms that 
underlie the TGF-P capacity in suppressing tumor development in 
normal cells, and conversely, to facilitate cancer progression and 
increase number of malignant cells. Our multi scale model focuses 
on the TGF-fi activation requirement, its autocrine properties of 
TGF-ji, its role in promoting cell contraction and being activated 
by cell contraction [6] . First we model the single cell and then the 
cell population by investigating the autocrine and paracrine of 
TGF-p in a generic epithelial tissue. We investigate the relation 
between the concentration of TGF-P and its receptors with respect 
to the stage of the cancer. In this perspective, we model the early 
invasion of few malignant cells in a healthy tissue and the different 
role of autocrine and paracrine TGF-fi secretion. The autocrine 
and paracrine behaviour are explored on the light of the 
evolutionary ecology principle of malignant and healthy cell 
competition. It is known that the bone tissue is the preferred niche 
of breast cancer colonization; we present a model of the role of 
TGF-P in bone invasion and alteration of bone tissue remodeling 
dynamics. 

In summary, here, we introduce a set of scale specific 
mathematical models that render the multi scale behaviours of 
TGF-P allowing us to describe the early breast cancer develop- 
ment and the initial condition of the metastasisation process by 



using a level of description familiar to biologists in order to 
encourage experimentation and hypothesis testing. 

Results 

The development and spread of cancer, although caused by 
driver mutations producing variations in gene expression and 
signaling disfunctions, involves cytoskeleton biomechanical chang- 
es that modulate cell dynamics at the tissue level. The 
development of a tissue mathematical model requires considering 
the TGF-P autocrine and paracrine behaviour of the cells. 
Therefore, we focused on the interface between intracellular and 
extracellular compartments. Given the different time dynamics of 
the reactions at the intracellular pathways level and the cell 
dynamics at tissue levels, we prefer to buUd distinct models, 
coupled by time averages of the fastest dynamics. Following the 
model developed by Laise et al. [14], we have focused on 
autocriny and paracriny behaviours of the TGF-P . Next, we have 
considered a tissue model to describe the effects of the TGF-P on 
cellular populations characterized by different driver mutations. 
Finally, we consider the bone niche model which allows us to 
describe the effects of the tumoral cells on the BMU (Basic 
Multicellular Unit) remodeling cycle. Each of these models 
describes different aspects of the TGF-P at a particular scale 
and they are loosely coupled by using averaged quantities of TGF- 
P in such a way to mimic the interactions between different scales; 
This allows us to consider each model as a "sub-model" which is 
part of a more comprehensive multi scale model. From here on, 
when we refer to the multi scale model, we will use the single term 
model, and we will specify otherwise; wh(;n referring to one of the 
scale specific model. Describing the dynamic development of 
tumors requires the knowledge of numerous degrees of freedom, 
which often are not experimentally available, and a complex 
model able to correctly analyse all the data, retrieve the relevant 
information of the present state and simulate its future evolution. 
Here, we propose a model to explain the early stages of tumor and 
its evolution in bone tissue based on production and sensing of 
TGF-P in both the paracrine and autocrine processes. 

The Intracellular Generation of TGF-fi Fluxes 

In this work, we follow the model of Laise et al. and [15,16] that 
for the sake of completeness, we will re-propose and adapt into the 
structure and the arms of a multi dimensional model that embrace 
both the intraceUular/ceU and cell/ tissue interfaces. They 
considered a simple biochemical model which explicitiy describe 
the key features of the Smad pathway. The first three equations of 
the model, see system of equations (1), accounts for the interaction 
between TGF-P {fi^ and the inactive membrane receptor (i^ec)- 
Following a successful encounter, the receptor turns into its active 
form here denoted 7?*^. In reality, different isoforms of TGF-P and 
of its receptor exist. In [14], the authors simplify the model setting 
by considering just one receptor type, which can operate in its 
active (bound to TGF-P) or inactive configuration. The activated 
receptor is able to phosphorylate the i^-Smad proteins (^c) in the 
cytoplasm, resulting in the formation of a new species, the 
phosphorylated Smad protein, here labelled p^c. Once phosphor- 
ylated, the Smad proteins head towards the nucleus. The 
translocation of the Smad proteins into the nucleus (p5n) is 
necessary to a( ti\ ate the transcriptional activity. This is a complex 
process, possibly organized in cascade regulatory cycles. Note that 
cytoplasmatic pS^ are modified into nuclear with a rate 
specified by the control parameter k\^. The presence of 
phosphorylated Smad in the nucleus is in turn associated to an 
increase of the hmgal gene expression [17]. The nucleus is also 
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Figure 1. 3D representation of the mammary duct. The mammary duct is formed by epithelial cells. Normal ductal cells (i^ = 0) are regularly 
arranged on a single 2-dimensional layer. Tissue irregularities appear in presence of benign ductal cells, but cells still lie on a surface. When epithelial 
cells develop an aggressive phenotype ((l>^3), they begin to stratify. 
doi:1 0.1 371 /journal.pone.0088533.g001 

enriched by phosphates (PPase) targeting phosphorylated Smad 
[16]. The dephosphorylatioii of the nuclear Smad proteins pSn 
results in non-phosphorylated nuclear Smad elements Sn- The 
model is composed of eight variables, respectively P^^, R^c, R^^, S^, 
pSc, pSn, S'n and hmgal. The associated concentrations obey to 
the following set of ordinary differential equations: 



5/i5ec=-'-b PeciORecit), 

8,R,c=-rbPJt)R^it), (1) 

d, i^:,=;-b ;8ec(0 Rec{t)-kp RUt) Scit), 



PLOS ONE I www.plosone.org 



3 



February 2014 | Volume 9 | Issue 2 | e88533 



Modeling TGF-p in Cancer Tissue Dynamics 



8, j,S, = kp R*Jt) 5,(0 --tin pSdO, 
8t pSn=k[^ pScit) — ks p>S'n(f)— ^dephos p'S'nCO ^^ase, 
8t Sn^kdsphoa pSnit) PPase-kex Snit), 
8i hmga2 = km pSn{t). 



where the notation 8t represents the partial derivative in respect 
to the dependent variable / indicating tlic \'ariati()n in time of the 
following quantity on the right. In the last equation of the system 
(1), hmgal measures the mRNA amount and as such assumes 
arbitrary units. All the other concentrations are instead expressed 
in nmol. Notice that the concentration of PPase remains 
unchanged all along the dynamical evolution. Asymptotically, 
the system admits a family of stable fixed points. From the 
experimental results in [16], eventually the TGF-fi amount as well 
as the activated receptors amount go to zero, and consequentiy, 
there is only one fixed point that satisfied such conditions. 
Analogously, according to the model specification, we have 
i^ec = ^ec = p'^c = p'^n = "Sn = 0. The quantities R^z, Sc and 
hmgal converge instead to stationary solutions, function of the 
initial condition and of the kinetic parameters involved. The TGF- 
P pathway is in particular reduced to a limited set of meaningful 
chemical reactions that are presumably implicated in the 
transmission of the signal from the cell surface, as triggered by 
TGF-fi, to the cell nucleus. The analysis conducted by the authors 
in [14] is in particular aimed at inspecting the out-of-equihbrium 
dynamics of the system, as driven by the externally imposed TGF- 
P. The model proposed by Laise et al. for the epithelial- 
mesenchymal transition predicts the concentration of mRNA 
associated to gene hmgal, properly describes the results of the in- 
vitro experiments set up by [16] and reproduces the right 
unperturbed steady state characterized by specific concentrations 
of cytoplasmatic Smad proteins and unbound receptors R^, 
which have been carefully evaluated by Schmierer and collabo- 
rators. 

To address the problem of building a model which takes into 
consideration the effects of the TGF-f} pathway signaling and 
tumor regulation at different scales, we adopt this intracellular 
model as a starting point. In our multi scale approach, we 
introduce two main re-adaptation of the previous intracellular 
model. First, we have done a model order reduction regarding the 
Smad and the hmgal concentrations by imposing kp constant 
in time, which allows us to simplify the intracellular processes 
disregarding the effect of the Smad signaling cascade reactions, see 
equations (1) in the dashed line box, in order to focus on the TGF- 
P pathway, see equations (1) in the continuous line box. The 
justification for such approximation, as explained in [14], is due to 
the large amount of Sc in comparison to R*^ and to the re- 
integration of Sc which permits us to consider the variation of the 
Smad in the cytoplasm negligible. Second, for a description of 
processes that occur not only at difiFerent spatial scales, but also at 
different time scales, we need to introduce source terms (and sink 
terms when necessary) for the synthesis of both the TGF-P and its 
receptors in order to move the fixed point in such a way the 
remaining quantities, at intracellular scale, are all different from 
zero. This is much more reasonable for long and continuously 
ongoing processes; while positive quantities, like the concentra- 



tions, approaching zero implies the processes will come to a 
standstill irremediably affecting the system at all scales. 

The Cellular Model 

The cellular model is the first of the three models where we 
discuss the intracellular/ cell interface and take into consideration 
the dynamic effects of the TGF-p in a small patch of cells. With 
this model, we address the problem of production and internal- 
ization of TGF-p along with the binding of autocrine or paracrine 
TGF-P to the receptors on the membranes without considering 
the detailed spatial disposition of cells. Even though the autocrine 
and paracrine signaling are completely distinct forms of exchang- 
ing chemicals, it is impossible to distinguish between the two when 
they occur at the same time. It is true that a cell can sense the local 
spatial inhomogeneity of chemicals and the heterogeneity of and 
positions of other cells it is in contact with; hence, the cell can 
regulate itself to secrete the chemical compounds along preferred 
directions in such a way the chemicals will most probably follow an 
autocrine pathway or a paracrine one. Nevertheless, these cellular 
behaviours and this level of detail are unknown and unavailable 
for the TGF-p signaling. 

On the other hand, it is important to stress that, if all the cells 
have similar behaviours in respect to the TGF-P secretion/ 
absorption, and they are homogeneously distributed, then the 
average paracrine signaling to and from other cells cancel out one 
another. In the case proposed here, the paracrine signaling inward 
and toward the neighbour cells do not cancel out because cells 
with different phenotypes do not behave in the same way. To 
overcome such difiiculty, we propose to subdivide the space where 
cells are placed by adopting a "cell-centric" point of view. 
Therefore, for each cell, we consider a spheric volume containing 
and centred on the cell itself. The surface of the cell membrane 
divides this volume between the intracellular part and the external 
part. The latter is divided in the extracellular region, which is the 
closest neighbourhood just around the surface of the cell 
membrane where the TGF-P produced by the cell is released 
and can bind with the receptors of the specific cell at the centre of 
the volume, and the diffusive region representing the farthest part 
of the spherical volume from the cell where the TGF-P cannot 
reach the receptors on the cell surface of any cell. Between two 
nearest cells, the farthest parts of the respective diffusive regions 
are the intangible frontiers where the two "cell-centric" spatial 
subdivisions join together. Furthermore, we extend the previous 
definition by considering the diffusive region of each cell as the 
place where all surround "cell-centric" spaces join together; 
therefore, the diffusive region is the common area between a cell 
and all its nearest neighbour cells. Tlu- TGF-P entering in the 
diffusive region loses the possibility to bind with cells and also loses 
any dependency on the cell from whom it has been produced. In 
other words, the TGF-P produced by the cells, which does not 
bind autocrinely, flows before into the diffusive region and then 
flows back toward all the cells in the neighbourhood indistinctiy, so 
that it can bind in a paracrine way. 

The portion of volume occupied by the union of all the external 
parts depends on the volume of the cells and the distances between 
their membranes. This intercellular space is filled with the 
extracellular matrix (ECM), a fibrous mesh which gives it the 
peculiar behaviours of the porous media. The ECM is responsible 
for the reduced diffusivity of the active TGF-p. As a consequence, 
the TGF-P in the diffusive region cannot easily diffuse toward the 
cells, and the TGF-P in the proximity of the cells is more easily 
conveyed toward the receptors on the cell surface. The 
prolongated time of diffusion is also responsible for the dispersion 
of TGF-p. The obstructions in the porous media introduce non- 
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linearity in the dififiision equation. Excessive non-porosity can 
forbid the diffusion of molecules through the media, or it will 
oblige molecules to foUow highly crooked paths. To avoid 

excessive complications in the diffusion of TGF-p due to the 
complex geometry produced by the disposititm of the fibrous 
mesh, we can take into consideration tlu- a\crag(-d conformational 
characteristics of the ECM in the diffusive region by using an 
effective diffusion coefficient which can be experimentally 
measured. Furthermore, the compartmental description of the 
ECM as a region with impaired diffusivity also reproduces the 
function of storing growing factors. 

Using the first Pick's law [18], the flux of molecules of TGF-fi 
crossing a unitary orthogonal surface is given by 
J= —Def[8x(P/V^[) where the effective diffusion coefficient D^g 
includes the porosity of the ECM, and the density of TGF-P is 
expressed as the quantity of molecules over the diffusive region 
volume Fdf . Discretizing the first Fick's law and multiplying both 
sides of the equation by the orthogonal surface ^df crossed by the 
TGF-P during the diffusion between the extracellular regi{m and 
the diffusive region, we get JA^ = — .Deff^df / f^df x A/V ^-^ where 
Viif / Aijf » Ax-lOvm is the distance between two adherent cells 
whose neighbouring membranes are bridged by E-cad. The 
quantity of molecules escaping through the orthogonal surface per 
unit of time can also be expressed in terms of the paracrine rate 
JA^f = — rp^Ap and hence, the TGF-P paracrine rate is equal to 

We want to remark that the introduction of the diffusive region 
makes it possible to mimic the paracriny between cells, but at the 
same time, we completely disregard the relative disposition 
between cells. Such approximation is not entirely appropriate to 
describe the exchange of TGF-P between two specific cells. 
However, it is justifiable when describing an averaged effective 
exchange of TGF-P between different groups/ types of cells. 

After the TGF-P binds with the receptors, a series of reactions 
involving the Smad proteins foUow their internalization. For the 
sake of simplicity, we do a model order reduction on the 
intracellular model proposed in [14] by considering the non- 
phosphorylated Smad proteins in the cytoplasm constant in time, 
see last equation in the continuous hne box and first equation in 
the dashed line box of the system of equations (1). 

Due to the importance of TGF-p for different aspects of the 
ceUular life cycle, degradation of TGF-fi in the intracellular 
compartment can not be neglected especially in healthy cells 
where the over-accumulation of TGF-fi can produce a large 
disregulation. Therefore, we consider that healthy cells ubiquiti- 
nate part of its TGF-p, while mutated ( (-lis do not perform such 
activity. The set of partial differential equations for the cellular 
model are: 



rbiSec(0>O-Rec(fO + '-pc 



/?df(0 



(2) 



d, R^i<l>,t) = r,y^ [RsM(<l>)-R^((t>,t)]-rt, P^(^,t) R^ilt), (3) 



a, R*^i<^,t) = rb P^ilt) R^(^,t) - [rsgn + ^^,o ru] R*J^J), (4) 



StPdf(t) = rpc 



(5) 



The .system of equations (2-5) defines the evolution in time of the 
TGF-p produced [P^^, its receptor on the cell membrane (Rcc), the 
internalised TGF-fi [R^^ and the TGF-fi in the diffusive region 
(jSdf)) see Table (1). While the TGF-fi in the diffusive region 
represents the total amount of TGF-fi that all the cells in the nearest 
neighbourhood are paracrinely exchanging, the other variables are 
intended as averages all over the sub-populations of cells with the 
same phenotype. Therefore, the quantities fi^, R^c and R^ depend 
on the phenotype expressed by the index (j> (and (p). The region of 
interest within the tissue where the tumor begins to develop is also 
partitioned in sub-populations of cells, p, identified by the 
phenotype index which goes from 0 to a maximum value 3>. The 
phenotype index is associated with the cell aggressiveness and 
sensibility to the TGF-fi so that (^ = 0 corresponds to the healthy 
cells and as (j) increases, the more cells are aggressive and need more 
TGF-fi in order to respond to its signal. 

In equation (2), the first two terms describe the synthesis of TGF-fi 
that is secreted in the external region where it is activated and the 
binding between the TGF-fi and its receptor. The third term takes 
into account the averaged values of TGF-fi transferred by one cell 
with phenotype ^ toward to the diffusive region and the mean TGF- 
fi per cell received from the diffusive region. Similarly, equation (3) 
describes the synthesis of TGF-fi receptors which are displaced on 
the cell membrane and binds with the TGF-fi present in the 
extracellular region. The TGF-fi binds to its receptor on the 
cytoplasm membrane, and it is internalized. Inside the cell, the TGF- 
fi interacts with the Smad [14] at rate r^gn, and to avoid an excessive 
abundance of this protein, ubiquitination occurs with rate r„, see 
equation (4). The operator Sjj in equation (4) is a delta of Kronecker 
which takc's value 1 when the indexes ;' = j and zero when the two 
indexes are different. The variati(m of total TGF-fi in the diffusive 
region is due to the incoming TGF-fi which each cell exchanges 
paracrinely and the out coming flux shared among all the nearest 
neighbour cells, first and second term in equation (5) respectively. 

We have used the TGF-fi and TGF-fi receptors gene expression 
data in the (x-Uular model equations (2-5) to evaluate the 
respective synthesized quantities and the concentration of the 
active TGF-fi in the intracellular compartment for each sub- 
population with phenotype index in the range [0,0]. In Figure (2), 
we show the results for different TGF-fi isoforms. Due to the 
disregulations of TGF-fi production in cells with driver mutations, 
we can observe that the concentration values of TGF-fi fluctuate 
between healthy, pre-tumoral and tumoral cc-llular phenotypes, 
and there is no clear trend which foUows the severity of the tumor. 
On the other hand, some particular cases are characterized by a 
monotonic relationship between the populations' phenotype and 
the concentration of the TGF-fi quantities, see the i?ec for the 
isoforms TGFB1I1/TGFBR3, both R*^ and fi^^ for the isoforms 
TGFB1/TGFBR2, the R* for the isoforms TGFBI/TGFBR2 
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and the Rec for the isoforms TGFB3/TGFBR3. These particular 
isoforms are usefitl to identify the tumor cell phenotype. For 
example, the expression for the TGFBIII receptors, which have 
an antitimiorigenic effect, decreases with the increase of cellular 
aggressiveness. 

It is important to stress that the dynamics of the TGF-ji pathway 
of a sub-population with a specific phenotype described by the 
system equations (2-5) depends on the dynamics of the TGF-fi of 
the other ceUular sub-populations in two different ways. The first 
one is the paracrine exchange of TGF-fi with rate Vpc between the 
sub-populations expressed in equation (2) and equation (5). The 
second way the equations (2-5) are coupled between different 
phenotypes with one another is given by the average number of 
nearest neighbour cells p{^,t). The former describes a cellular scale 
phenomenon, while the latter is a tissue scale phenomenon. Hence, 
if the number of cells of the different sub-populations is not constant, 
but changes dynamically, then it is necessary to supply the cellular 
model with a set of equations for the cellular evolution of the tissue 
(see next subsection for the tissue model). 

The Tissue Model 

Mutations are responsible for the behavioural changes of cells 
that, from a healthy state in which they are capable of correctiy 
sensing and responding to the surrounding signaling, enter to a 
mutated state, where the cells cannot self-regulate in response to 
the homeostatic signals. Mutations induced by external agents, or 
due to the occurrence of casual variations in the DNA' s 
transcription while proliferating, can be easily accumulated during 
a cell life or in multiple progenies. Hence, the increasing in the 
population number and the survival of mutated cells can obstruct 
the tissue integrity and its functional activity well before the cells 
acquire a highly malignant phenotype. As previously stated, the 
TGF-p signaling has a pivotal role in maintaining homeostasis at 
the ceUular scale and the functional integrity at the tissue scale. 



Table 1. List of variables. 




variable/unknown 


symbol 


initial value 


extracellular TGF-ll 


Ac 


{2,2,2,2} [nmoi] 


receptors on the 
membrane 




{2,2,2,2} [nmol] 


internalized receptor- 
ligand 


K 


{1,1,1,1} [nmol] 


autocrine TGF-ji 




0.1 [nmol Fjf-'] 


number of nearest 
neighbor cells with 
phenotype <j) 


P(<I>J) 


{6,0,0,0} 


osteocytes 


Ocy 


900 


osteoclasts 


Oc 




lining cells 


Lng 


0 [Sf„,-'] 


osteoblasts 


Ob 


0 [«f,ac"'] 


RANKL 


RANKL 


0 [nmol V,,,-'] 


BMP 


BMP 


0 [nmol F,.,^-'] 


CSF 


CSF 


0 [nmol Ff„,-'] 


BMD 


BMD 


100 [%| 


If the variable depends on the phenotype <j>, then we give a list of values sorted 
by increasing order of the phenotype. References and choices for the numerical 
values of the initial values are discussed in the data analysis section. 
doi:l 0.1 371/journal.pone.0088533.t001 



Indeed, the TGF-ji downregulates the cell proliferation, is 
responsible for the cell cohesion at high concentration, induces 
cell apoptosis [7]. Nevertheless, the anti-tumorigenic mechanisms 
provided by the TGF-fS, and their effectiveness, depend on the 
capability of the cells to properly sense its signal. DNA mutations 
(in particular, driver mutations [19]), on one side, can induce 
changes in cell phenotype which destabilize the correct cell 
functions. The weakening of a cell and loss of its stability are 
responsible for the increase of active TGF-p. On the other side, 
mutations can produce a resistance of the cell response to the 
TGF-P signaling. Cells with these driver mutations can also 
acquire the capability to produce a higher concentration of TGF-P 
which is required to reach a different homeostatic level without 
incurring in apoptosis. Eventually, a cell can always undergo a 
mutation resulting in the failure of the anti-tumorigenicity of the 
TGF-P and in an inversion of its role, meaning the transformation 
of TGF-P from Dr Jekyll to Mr Hyde occurs [8]. In the latter case, 
the TGF-P fails in downregulating cell proliferation and inducing 
cell apoptosis; while the excessive production of TGF-P becomes 
dangerous for the surrounding cells which have not yet acquired 
sufficient resistance to the growing factor. 

As previously said, cell mutations are random, and each 
mutation can induce apoptotic cell resistance to certain signals, or 
it can introduce ceUular instability and put the cell to death. 
Nevertheless, this does not mean there is no relation between 
mutations, or that a mutated cell can accumulate mutations and 
return to its original state. Indeed, the cell behavioural changes 
induced by mutations can be associated to a stage. This cell stage 
indicating the results of mutations does not variate continuously, 
but, on the contrary, it is characterized by abrupt phenotype 
changes. Therefore, mutations make the cell change from one 
stage to another with given probabilities, and the presence of 
correlations or constraints between mutations oblige the cell each 
time to variate its stage going through a small subset of all the 
possible mutation states. Indeed, the mutation state of a cell 
monotonicaUy increases during its life, producing a progressive 
change of its phenotype which from a normal stage goes through a 
series of pre-neoplastic steps to a neoplastic phenotypical stage. 

On the other hand, not all the occurrences of a mutation imply a 
change of the cellular phenotype or a change in the production of 
and response to the TGF-P . Furthermore, at each stage, a cell has a 
given probability in acquiring a complete resistance to the tumoral 
suppressor action of the TGF-P by switching to a phenotype where 
the TGF-P becomes a tumoral promoter. To address the diflFerences 
between cell tumoral behaviours and cell response to TGF-P, we 
introduce a discrete positive variable (j> that represents the cell 
phenotype. The index (j> subdivides cells into groups which share the 
same phenotype without considering the specific mutations 
accumulated by each single cell; therefore, each cell in a group 
has the same sensibility to the signaling induced by the active TGF- 
P bound to the membrane's receptors and activates the same 
amount of TGF-p per unit of time. From a biological point of view, 
in our prospective, the values of (j> has some relativeness with the 
hallmarks of cancer [2,3], and from a medical point of view, it may 
be related with the diagnostic stage of cancer. 

All quantities and characteristics which refer to a particular cell 
stage, or mutated cell population, wiU depend on the phenotype (j>. 
Therefore, normal cells have a phenotype </> = 0, pre-neoplastic cells 
correspond to = 1 , tumoral cells are indexed asfj> = 2, and cells with 
aggressive tumoral behaviours and strong resistance to TGF-P 
inhibiting signaling have phenotype 0 = cD = 3, [20,2 1] . A possibility 
of explaining the different responses to TGF-P signaling by cells 
with different phenotypes is thinking that a normal healthy cell has 
^ different mechanisms working in series to regulate the response/ 
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Figure 2. Table of TGF-ji synthesized, TGF-fi receptor syntliesized and TGF-P internalized for each phenotype c^. These values are 
obtained from the numerical solution of the cellular model equations (2-5) at time r=10 days. The various types of TGF-f} and receptors are 
expressed on top of each box. The tumor grades are sorted into increasing order of severity and aggressivity. 
doi:1 0.1 371/journal.pone.0088533.g002 
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sensing of TGF-fi. Therefore, at the beginning, the signal 
entering the cell is amplified to (y^Jj.)" where 0 < a < 1 and y is a 
suitable scaling constant. Then, the amplified signal is sent to all the 
Q> mechanisms and a binary information .se{0,l}, whicli indii atc's if 
there is an absence or presence of entering TGF-fi respec:tively, is 
sent to the first regulating mechanisms in the series. Each 
mechanism uses the entering signal to amplify the binary 
information 5 of a factor iyR^)". After all the <1> mechanisms are 
applied, the outcome (yJ^ej.)'*" compared with the amplified signal 
iV^tcT their ratio is used as an upregulation of the apoptotic 
signaling and as a downregulation of the cell prohferation. For a 
normal cell {(f) = 0), all the 3> mechanisms regulating the response of 
the TGF-jl are functioning, while for the successive phenotype 
^ = 1, there is one mechanism which always fails to function. The 
failure of the mechanism results in no contribution to the 
amplification of the binary signal. As the phenotype ^ increases so 
does the number of failing mechanisms. When all the mechanisms 
which amplify the entering signal are non-functioning due to the 
severe mutations like in aggressive tumoral cells, the binary signal s 
(non amplified) compared with the entering signal (jR*^)"' > 1 
produces a reduction/dumping of the TGF-fi signaling. Conse- 
quentiy, the cells respond by increasing the proliferation and 
decreasing the apoptosis. 

To describe the dynamic evolution of the cell populations for 
each phenotype (f>, we propose a model based on the effectiveness 
of the cell response to the TGF-P signaling. To easily describe the 
exchange of TGF-fi between cells and their phenotypical 
evolution, we focus our attention on a small region of healthy 
tissue in which cancer cells will form. Precisely, we consider a 
volume containing a cell, all its neighbour cells and part of the 
empty space inside the mammary duct, see Figure (1). Cells in a 
normal mammary duct are regularly arranged in such a way to 
form a tubular surface one cell thick. Even when cells start to 
present some pre-neoplastic phenotypical behaviours, they con- 
tinue to lay on the mono layer duct surface, but with a less regular 
arrang(;meiit. Therefore, after having sc't the \ ()lume of the region 
of interest, we can consequentiy fix the average maximum number 
of cells, C^ = Co = Ci = C2, which can lay on the surface of the 
mammary duct. On the other hand, neoplastic epithelial cells in 
the mammary duct tend to form multi-layers; therefore, the 
volume capacity, Gt, of this type of cells is bigger than in the other 
cases. The average number of cells p{ij>,t) within the volume of 
interest is described by the following equation: 

8,p((j>,t)=rp{(j>)\l -= —J 

-r,[YRU<l'Jr^*^]p(M-'-<iP{M 

f ®-i 1(6) 
+ rm< &^,<s>pis> X] Pif't) + (1 [(1 - '5^,o)( 1 -P<s>)p{^ - 1 ,0 - P(^,t)] I 



+ (1 -i5fo — '5f(^)- 



'•p(^-l)/', Y:%lp(.f>t) 



-U) 



-1)' 



with the additional condition that p((j>,t) = 0 for <j><0 (the 
phenotypes (fxO are set constantly equal to zero, because they 
do not have any biological meaning and they can be disregarded). 
The first term on the right hand side of equation (6) describes the 
proliferation limited by the volume capacity C^, and downregu- 
lated by both the TGF-fi entering the cell, R^, and the capabifity 



to respond to it, gpi(j>). The second and third terms express the 
cellular death induced by the TGF-fi signaling, which also 
depends on the phenotype sensing exponent gii(<l>), and by the cell 
instabilities induced by random mutations. The fourth term 
describes the changes of phenotypes (the increase of ^) as a 
consequence of the mutations. The delta of Kronecker shows that 
normal cells can only develop anomalous behaviours, and 
aggressive cells {(j> = <!>) do not change their phenotype. Similarly, 
the fifth term describes the occurrence of important mutations and 
of phenotype changes during the cell proliferation. We introduce a 
upper limit <J> for the phenotype values because when cells 
accumulate to many mutations, they reach a stage of instability 
which are inconsistent with both the aggressiveness of the cell 
phenotype and the diagnostic stage of cancer. A similar upper limit 
used to label the tumoral stage of the cancer cells has been adopted 
in [22] as the limit in which cells are prevalently characterized by 
an apoptotic regime instead of those characteristic hallmarks 
associated to cancer development [2,3]. 

In Figure (3), we show the numerical solution of the tissue model 
coupled with the cellular model. The figure represents the 
evolution of the average nearest neighbour cell sub-populations' 
densities inside the mammary duct. The occurrence of driver 
mutations inducing phenotype variations and the following 
disregulation of TGF-fi cell production results in an increase of 
cell populations with higher phenotype indexes to the detriment of 
the healthy cell population. Figure (3.A). In Figure (3.B) and 
Figure (3.D), we can see that pre-tumoral cells, having de\'eloped 
only a partial resistance to TGF-fl signaling, proliferate much 
slower than cells with aggressive tumoral behaviours. Nevertheless 
due to rare occurrences of driver mutations inducing phenotype 
changes, the cell population with phenot)'pe (j) equal to 2 increases 
much slower than pre-tumoral cells. Figure (3.C). In Figure (3.D), 
the tumoral cell sub-population with phenotype (j) equal to 3 is 
characterized by a long initial part where they are almost zero 
followed by a high proliferation phase until the maximum cell 
capacity is reached. This shows that, even though TGF-fi has lost 
its anti-tumorigenic role on aggressive cells by increasing their 
proliferation and driver mutations immediately changing the cell 
phenotype to the maximum value $ can occur at any moment, the 
delay with which aggressive tumoral cells form points out that they 
are prevalently originated by sub-populations which have already 
developed TGF-fi resistance more than by healthy cells. Hence 
this highlight the strong capability of TGF-fi in slowing down the 
cancer development by acting on pre-tumoral cells. On the other 
hand the steep increase of aggressive cancer cells, after the first 
plateau-like phase, remarks the role of TGF-fi as cancer promoter 
on aggressive tumoral cells. 

It is important to point out that both the cellular model and the 
tissue model describe the dynamic of averaged quantities. Also the 
coupling between the two is regulated by the a\'cragc of the TGF- 
fi entering the epithelial cells and the av('ragc numbers of 
neighl)our cells with a specific phenotyp(; which compose the 
mammary duct tissue in the region of interest. The differences in 
the order of magnitude between the rates in the cellular model 
with respect to tissue model, see Table (2), allow us to overcome 
such difficulties; in fact, phenomena occurring at the cellular scale 
are much faster and relax to a steady state more rapidly than the 
events happening at the tissue scale. Hence the difference between 
the two time scales and the loosely coupfing, due to averaged 
TGF-fi quantities, makes it possible to reduce the coupling 
between the extracellular and the tissue scales resulting in multi 
scale model consistent with the biological phenomena. 

Different factors play a role in the prediction of cancer 
evolution. The presence of characteristic time scales in the 
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A 

5.05 r 




Figure 3. Evolution in time of tlie average sub-populations' 
densities. Given a volume Vnagh of a size equal to the average volume 
that, in a healthy mammary duct tissue, contains exactly a cell and all its 



nearest neighbour cells arranged on one layer, the densities for each 
sub-population In the volume Kneigh is shown. Starting from the top. A) 
the healthy cells, B) the pre-neoplastic cells, C) the tumoral cells and, D) 
the aggressive tumoral cells are graphed respectively. On the abscissa 
the time in months. 
doi:l 0.1 371/journal.pone.0088533.g003 

development and the dynamics of breast cancer is not so easy to 
define or to detect. Nevertheless, few considerations can be made. 
Observing the values of the parameters of the multi scale model, 
some of which are given in the literature and others are chosen so 
that their order of magnitude is consistent with the range of values 
present in similar biological situations (see Sec. Methods), one can 
a posteriori pinpoint that there are two distinct time scales. The 
first describes processes occurring at intra-extra celltdar time scale 
and the other, at tissue time scale. Hence, the differences between 
the two time scales reflect the velocities at which the processes 
happen at both the cellular level and the tissue level. This implies 
that fluctuations due to various cells activities, e\ en when they are 
large, rapidly decay and cancel out. From the tissue scale point of 
view, we can disregard large rapid cell variations, and conse- 
quently, a posteriori, we are justified in adopting the mean-field 
form of interaction between the cell and the tissue scales. 

Looking at various clinical cases and at their respective gene 
expressions, these two time scales cannot be explicitly determined 
due to the lack or impossibility of direct experimental measure- 
ments. On one hand, using the information for each patient to 
adjust the model parameters so to better fit each dynamics of the 
disease, we can observe that there are two distinct time scales, even 
though they are not unequivocally defined. On the other hand, 
these two time scales do not identify different steps of the cancer 
development, nor do they represent two different dynamics. They 
only express that an abrupt event which suddenly changes the 
behaviours in a small group of cells may initiate the disease, but its 
effect will be felt at larger scales only at longer times thanks to the 
TGF-p which represents the vector of the interaction. 

The bone Niche Model and the Effect of TGF-p In Bone 
Remodeling 

Breast cancer bone metastases are predominantly osteolytic and 
accompanied by bone destruction, bone fractures, pain, and 
hypercalcemia, causing severe morbidity (bone metastases occur in 
about 70% of patients with advanced breast cancer). Comorbidity 
addresses the occurrence of different medical conditions or 
diseases, usually complex and often chronic ones, in the same 
patient. There is a causal effect, and the bone is a unique 
microenvironment in which breast cancer thrives [23]. Bone is 
continuously being formed by osteoblasts and resorbed by 
osteoclasts, not only to maintain mineral homeostasis but also to 
cope with the microfractures that occur naturally [24—27]. We 
belie\e that txjmputational modeling could be very effective in 
shedding a light across the intrinsic difficulties of integrating 
evidence obtained from experiments and observations spanning 
different scales of time and space. There is a growing number of 
mathematical and computational models investigating the com- 
plexity of this dynamics [28-,34] and the interaction with cancer 
cells [35,36]. In the adult skeleton, TGF-fl is abundant in the bone 
matrix, where is released following the initiation of resorption 
TGF-fi is released from bone matrix [6,37,38]. Few recent studies 
have highlighted the complexity of breast cancer metastasis in 
bone microenvironment [6,9,39,40]. Although TGF-fS enhances 
the recruitment and proliferation of osteoblast progenitors, TGF-fi 
potently inhibits later phases of osteoblast differentiation and 
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maturation and suppresses matrix mineralization [41^5]. Oste- 
oblasts are derived from mesenchymal stem cells and their primary 
function is to synthesize the organic collagenous matrix and 

orchestrate its mineralization by producing bone matrix proteins 
including osteocalcin, osteopontin and bone sialoprotein, and 
providing optimal environmental conditions for crystal formation. 
FuUy diflFerentiated osteoblasts that are surrounded by mineralized 
bone tissue are called osteocytes and act as mechanosensors in 
bone tissue. They are the most numerous cells within the bone 
tissue and scattered evenly through the matrix. With their 
flattened morphology and long processes, they form a sensory 
network which allows the detection of abnormal strain situations 
such as generated by microcracks. Osteocytes are connected to 
one another and to surface osteoblasts (bone lining cells) via gap 
junctions. Osteocytes secrete sclerostin which is a master switch to 
prevent the body from making too much bone. The tuning of 
sclerostin allows osteocytes to control the osteoblast's activity in 
bone formation. Under physiological mechanical stimuli osteo- 
cytes prevent osteoclast's bone resorption by changing the 
RANKL/osteoprotegerin (OPG) ratio. By communicating these 
signals to bone lining cells (the second terminally differentiated 
osteoblast cell type) or secrete factors that recruit osteoclasts, 
osteocytes initiate the repair of damaged bone. We described a 
model of the interplay between osteocytes, osteoblasts and 
osteoclasts. This model focuses on the release of TGFP during 
resorption phases. Pathological conditions can alter the equilibri- 
um between bone resorption and bone formation. Often, 
disregulations fa\'ouriiig osteoclastogenesis; osteoporosis is an 
example of negative remodelling: the resorption process prevails 
on the formation one and this reduces bone density, so increasing 
the risk of spontaneous fractures. Here we model the release of the 
TGF-P from the bone matrix upon the action of the osteoclasts 
which favour the breast cancer metastasis and the unbalance of the 
bone dynamics. While experimental works represent primary 
sources of parameter values [46], the mathematical and compu- 
tational recent works, such as [28-30,32], provide a valuable 
validation and discussion of the range of parameters value. Here 
we have used parameters values accepted by various Uterature. 
There is a growing use of omics genome wide analysis, see [47] 
among others. In our work we have estimated some parameters 
from gene expression data. It has been pointed that the difference 
between an approximate and exact model is usually remarkably 
smaller than the difierence between the exact model and the real 
biological process [48] . Taking into account the recent experimental 
results on the central role of osteocytes [24] and the most recent 
mathematical models, in particular Pivonka [43] we have developed 
the following conc(;ptual model: first we stress the importance of the 
osteoblasts forming the lining of cells attached to the bone tissue. It is 
likely they are in communication with the osteocytes in the bones via 
a network of canaliculi. When osteocytes die because of fracture 
there is a loss of communication of the lining osteoblasts with the 
osteocytes. This happens together with a release of RANKL whit:h 
is a signal for osteoclasts to arrive and it is a diflferentiation signail for 
the osteoblast which may take sometime to happen. We have 
modeled this delay by using a delay differential equation. The 
osteoclasts will carry on a resorption of the bone which could be 
stopped by decrease of RANKL or by the release of TGF-jj by the 
bone. Then the osteoblasts with start deploying bone material and 
proteins and complete the differentiation process by rebuilding the 
network of canaUculum processes. One can think that cancer cells 
would take the advantage of absorbing part of the TGF-fJ so 
decreasing the apoptotic probability of the osteoclasts. We have 
used DDE to implement both the differentiation lapse of the 
osteoblasts. Delays not only are an explicit representation of the 
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time necessary for spatial transfer of information, but are important 
also to take into account the presence of other reactions which are 
not exphcitiy expressed, or included and the lapse of time for them 
to happen. The osteoblasts during the bone resorption and bone 
formation appears in three forms: uncommitted osteoblasts, pre- 
osteoblasts and active osteoblasts. The differentiation from the first 
type and the second is accelerated by the presence of TGF-fi and 
the differentiation from the second type to the third is inhibited by 
the TGF-fi. Uncommitted osteoblasts lay on the surface of the bone 
vessel and their differentiation process is influenced by the 
osteoclasts' resorption activity which is mainly triggered by bone 
fracture (and microfracture). This event propagates and after an 
average time To uncommitted osteoblasts begin to differentiate. The 
presence of TGF-p can be of help for their differentiation, or even 
necessar)-, but the TGF-fi cannot abbreviate it because it does not 
speed up the travel of information. On the other hand TGF-ji can 
delay the activation of osteoblasts. The complete and complex 
differentiation process of the osteoblasts can be summarized from 
the first stage to the last stage by packing all the complexity of the 
process inside a delay ttotiP) = To -l- T^diffifi) dependent on the TGF- 
fi, p. The delay idiff(i^) must be positive and finite, because the 
TGF-fi does not preclude the osteoblast activation and second 
important factor to remark is that the effect of TGF-fi on the 
osteoblasts is local (at contact) in space and picked (rapid effect and 
defined delay ) in time. On the other hand the diffusion of TGF-fi, 
osteoblasts and their concentration are much more sensible to space 
inhomogeneities and large time distributed events. Our DDE model 
equations are listed below: 

d, Ocy =focy Ob{t) - aocy ©(f) 0(1 - 1\ (7) 
Oc = roc RANKL{t) CSF{t) - oqc Oc{i), (8) 



d, Lng=pLng BMP{i)-mL„g BMP{t-x{P)) 
-aLng Lng{i), 



d, Ob=mLns BMP{t-x{Pj)-vfocy Ob{t)-aob Obit), (10) 



d, RANKL=prkl {Ocy-Ocy{i)) ®{t) 0(1 -f) 
+rRKL Ob(t)-dRKL RANKL, 



d, BMP=pbmp {Ocy-Ocy{t)) -dsup BMP(i), (12) 
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Table 2. List of parameters. 





parameter 


symbol 


value 


gene expression of the TGF-/^ 


PsM 


{-1.853, 1.724, 1.614, 1.669} [nmol] 


gene expression of the receptors 




{-7.47, 6.713, 5.017, 4.128} [nmol] 


rate of synthesis 




133.333 [s"'] 


Smad signaling rate 




0.003018 [s"'] 


paracrine rate 




10 [s"'] 


binding rate 


fb 


0.0074 [nmor' s"'] 


ubiquitination rate 




0.007 [s"'] 


total number of phenotypes 
and malignant phenotype 


0 


3 


proliferation rate 




{-1. xlO""*, 1. xlO"^ 1. xlO"^ 1. xlO""^} [s"'] 


apoptosis rate 




9. xlO"" [s"'] 


degradation rate 




1. xlO"^ [s"'] 


mutation rate 


I'm 


1. xlO"" [s"'] 


proliferation exponent 




{-0.12, 0.06, 0.,-0.06} 


apoptosis exponent 


n (S) 


{-0.12, 0.06 0. —0.06} 


probability of malignant mutation 


Po 


5. xlO"' 


cell capacities 


Co 


{-6, 6, 6, 8} [Vneigh"'] 


dimensional constant 


a 


1 [5-^1 


dimensional constant 


Y 


1 [nmol"'] 


Osteocytes formation rate 


focy 


0.00032 [Sfrac l/osteon"' day"'] 


Osteocytes apoptosis rate aOcy 


^Ocy 


30 [l/osteon"' day"'] 


Osteoclasts recruitment rate 


'Oc 


0.8 [Vfrac nmol"^ day"'] 


Osteoclasts apoptosis rate 


^Oc 


0.0625 [day"'] 


Lining cells proliferation rate 


Ping 


50 [jim nmol"' day"'] 


Lining cells maturation rate 




50 [^/m nmol"' day"'] 


Lining cells apoptosis rate 


^Lng 


0.1525 [day"'] 


Osteoblasts apoptosis rate 


<^Ob 


0.1525 [day"'] 


RANKL production (Osteocytes) 


Prkl 


300 [nmol Vosteo^ Vfr^r' day"'] 


RANKL production (osteoblasts) 


^RKL 


2.5 [nmol S,,,^ l/fra^"' day"'] 


RANKL degradation rate 


^RKL 


8.6643 [day"'] 


BMP production rate 


Pbmp 


1 [nmol Vosteon Wrac ' day"'] 


BMP degradation rate 




8.6643 [day"'] 


CSF production rate 


PCSF 


0.1 [nmol l/„s,eon V„ac"' day"'] 


CSF degradation rate 


^CSF 


69.3147 [day"'] 


Bone resorption rate 


'bn 


9.4x10"" [l/f,ac nmol"' day"'] 


Bone formation rate 


fbn 


5.9098x10"*^ [Sfrac nmol"' day"'] 


TGf-fi extraction rate 


'K 


2.2x10"'° [nmol day"'] 


TGF-/i degradation rate 


df 


0.1189 [day"'] 


Number of cancer cells 


Nc 


2000 


Constant time delay 


To 


6.02 [day] 


Constant time delay 


A 


22.36 [day] 


Maximum amount of osteocytes Ocy 




900 [l/os,e„„"'] 


Minimum effective amount of TGF-/i inducing delay 


£ 


8x10"'' [nmol 1/frac"'] 


Burried osteoblast factor 


V 


1 [^osteon Wrac ^] 


If the parameter depends on the phenotype then we give a list of values sorted by increasing ord 
values of the parameters are discussed In the data analysis section. 


2r of the phenotype. References and choices for the numerical 



doi:10.1371/journal.pone.0088533.t002 
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e, CSF=pcsF {Ocy-Ocy{t)) -dcsF CSF{t), (13) 



d, BMD =fbn Ob{f) - Tbn Oc{f), (14) 



= rjs Oc{f) - dfs m - N, V m, ( 1 5) 



The equation (7) describes that a microfracture occurring at 
t = 0, is sensed by osteocytes that undergo aj)optosis \\-hi<:h ends in 
about one day [46,49] . In equation (8) the osteoclasts are recruited 
to the BMU in response to a combination of RANKL and CSF, 
and die at a rate Uqc', The equation (9) shows that immature 
osteoblasts are recruited in response to TGF-fi here represented as 
BMP, and differentiate into mature osteoblasts after 20 days. 
Mature osteoblasts can either self-bury with a rate vfocy and 
differentiate into osteocytes at rate focy or die at rate- aob, see 
equation (10). References could be found in [46]. The equation 
(11) describes that surviving osteocytes secrete RANKL at a rate 
proportional to the "size of the fracture" i.e. number of osteocytes 
that underwent apoptosis; osteoblasts produce both RANKL and 
OPG at rates dependent on their maturity and with characteristic 
delay [36]. The equations (12, 13) describe that BMP and CSF are 
secreted by osteocytes in response to the microfracture damage. 
Both chemicals decay, but CSF is also depleted by osteoclasts; 
representing the idea that CSF is a chemoattractant and the 
osteoclasts move on when its concentration becomes low [50,51]. 
Note that TGF-fi and bone morphogenetic protein (BMP) have 
often opposite behaviour. In a hair follicle stem cell niche BMP 
signaling maintains stem cell quiescence while TGF-P stimulates 
stem cells both in vivo and in vitro and antagonizes repressive 
BMP signaling [52]. 

The equation (14) follows [28] equation where bone is rebuilt at 
rate fi,„ by mature osteoblasts and resorbed at rate ri,„ by 
osteoclasts. The equation (15) describes the TGF-fi from the 
resorption action of the osteoclasts. 

In Figure (4) we show the simulated results of the ODE model: 
Bone density (A, B), number of osteoblasts (C, D), number of 
osteoclasts (E, F) and TGF-ji (G, H) compared between control (A, 
C, E, G) and cancer condition (B, D, F, H). Results show a 
negative remodelling balance in the cancer case. The values in 
Figure (4) are expressed in terms of averaged quantities which are 
the osteon volume, the fracture volume and the fracture surface. 
Considering a microfracture of 1/im width which caused the 
damage of a number of osteocytes 30 times smaller than the 
number of osteocytes in the osteon, the relation between them are 
^osteon ^BOFfrac and Ffrac-^ftac-l/^m, see [5,3]. 

The negative balance of the bone density matrix due to the 
effect of cancer cells in the bone niche, where metastasis occurs, is 



strongly depending on the number of cancer cells. Even though 
breast cancer cells find in the bone tissue a richer environment of 
TGF-fi favouring their proliferation, the formation of metastasis is 
not a very easy and probable event. Different biological defensive 
systems and causes concur to avoid the formation of metastasis. 
Nevertheless maybe due to the high number of breast cancer cells 
detaching from the main tumor and reaching the bone tissue 
though the vascular system, or maybe due to the occurrence of the 
rare event in which all the metastatic defensive systems fail, few 
cancer cells (not necessarily in positions close one another - 
information neglected in the models proposed) initiate the 
metastasis. In the very initial development stage, the few metastatic 
cells do not really affect the bone system. On the other hand, at 
this stage, a microfracture and the successive remodeling process 
proceed faster than the increasing of the local cancer cell 
population. In the early stages of metastasis formation process, 
we assume that during the remodeling process, the number of the 
tumoral cells in the bone niche is constant, but after successive 
microfracture, due to the abundance of TGF-fi, the cancer cell 
population increases. The more cancer cells increase the more 
they affect the bone remodeling process. When cancer cells 
increase in number, on one side, they increase the time of 
maturation of lining cells into osteoblasts through the action of 
TGF-fi on the BMP, on the other side, they obstruct the tissue re- 
mineralization. The result is a weaker and less dense bone. 

he TGF-fi is responsible for the delay in the maturation of 
osteoblasts. An excessive quantity of TGF-fi released during the 
bone remodfJing process is the cause of the reduction of the bone 
density. In the same way, a reduced quantity of TGF-fi induces a 
rapid maturation of osteoblasts. This prevents osteoclats resorption 
in the BMU causing a local increase of bone mass with less 
structural strength. 

Absence of a spatial representation of the bone niche does not 
allow us to completely describe the dishomogeneities in the bone 
mineral density caused by the cancer cells in bone tissue which are 
known as mixed lesions. Cell-to-cell spatial interactions like 
volume exclusion and chemotaxis are necessary to reproduce 
mixed lesions. In order to mimic the presence of mixed lesions, we 
can use our model to simulate the occurrence of multiple fractures 
in which the intensity of the released TGF-fi fluctuates randomly. 
The variability of the bone density accumulated at the end of each 
remodeling by spatially independent BMUs can be considered as 
an index for mixed lesions. Hence, high variability will be an 
indication of mixed lesions, while low variability will represent the 
cases of osteolytic or osteoblastic lesions depending on negative or 
positive changes of the averaged bone density, respectively. 
Nevertheless, such variability is related to the probability 
distribution used for the intensity of the released TGF-fi. 

In Figure (4), we show the evolution of the bone remodeling 
process after a small fracture when the population of cancer cells 
in the bone niche is of the order of 10^. This is sufficient to see that 
in the local region of bone in which metastasis is developing, the 
bone density matrix decreases of 1%. The bone mineral density 
represents a marker to identify the development of metastasis and 
its effectiveness is based on the causal relation or on comorbidity 
between breast cancer and osteoporosis. 

Discussion 

We have developed a mathematical model describing the multi 
scalar behaviours of the TGF-fi, first by representing the actions of 
the behaviours of the TGF-fi from the intracellular to the 
extracellular scale using a cellular model which is based on the 
results of the work pubhshed by Laise et al. for epithelial- 
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Figure 4. Bone remodelling process. Evolution in time of the bone density (A,B), the osteocytes (CD), the osteoclasts (E,F), the osteoblasts (G,H), 
and the TGF-ji In the bone niche released during the resorption (l,J). Control (first column) versus cancer condition (second column) is shown. The 
bone densities are expressed in percentage, the osteocytes are in number of cells per osteon volume (Kosteon), the osteoblasts are in number of cells 
per fracture surface (Sfrac), the osteoclast are in number of cells per fracture volume and the TGF-ji in nanomoles per fracture volume. On the 
abscissa the time is in days. 
doi:1 0.1 371/journal.pone.0088533.g004 



mesenchymal transformation [14] and which also addresses the 
important aspect of autocrine/paracrine signaling induced by 
the protein. Then we have further developed the model to 
investigate the extracellular mechanisms of action of TGF-p at 
larger scale and its involvement in cancer. This has led to a 
tissue model of breast cancer, which usually metastasise to the 
bone tissue. Therefore, we have also studied the relation between 
the metastatic cell, the bone tissue and the TGF-p. We have 
identified the TGF-P as one key molecule involved in cancer 
because it is antioncogenic and pro-apoptotic at low concentration 
and pro-tumorigenic at high concentration due to its mechan- 
obiology properties through the interaction with the cytoskeleton 
and with the extracellular integrins/E-cad proteins. This multi 
scale model is based on a multi scale multifunction molecule, in 
other words it is a sort of hub for apoptosis signaling at one scale 
and for the cytoskeleton function and cell-to-cell interaction at 
another scale. 

Our work has focused in generating and connecting three 
models: one cellular model (intracellular and extracellular TGF-fS) 
and two tissue models (the breast cancer in the mammary duct and 
the bone tissue). Each of these models refers to a specific level of 
organization, and hence it is possible to approximately associate 
each of them to a specific length scale. Such spatial subdivision of 
the models is also reinforced by their temporal subdivision in the 
dynamic processes (cell/tissue scales) and in the cancer disease 
progress (breast/bone tissues). The differences in temporal scales 
of the various processes described in these models hamper the 
aggregation of all these models into a single one. Indeed, the time 
scales of TGF-fi inside the cell are usually much shorter than 
events at the tissue level. This has been resolved by considering 
that the variables at the cellular level are passed as "averages" to 
the tissue level. The information related to the spatial positioning 
of cells in the tissue has been neglected in order to find results not 
dependent on the tissue geometry. In other words, we do not 
consider ceUular specialisation in the tissue; therefore, cellular 
parameters are space and cell type independent, but, if necessary, 
only phenotype dependent. In this way we avoid an increase in 
the number of degrees of freedom; furthermore, our approach 
leads to a model order reduction with no bias for the position of 
the cell. 

The TGF-j) has been reported to show an autocrine control 
system to inhibit cell replication, thus maintaining the tissue 
cell homeostasis. We have modeled the process that, following 
mutations, drives the inhibition into cell cycle progression. So 
the cell increases the autocriny i.e. the production and endocytosis 
of TGF-j) in order to regulate the ceU cycle progression, but this 
results in an increase of tumorigenicity; meanwhile, the TGF-j) 
also acts in disrupting the actin cytoskeleton which disrupts 
the actin and E-cad anchorage; hence, the cell contracts and, 
as consequence, this increases the activation of extracellular TGF- 
j); therefore the cell without contact inhibition starts replicating 
actively. The end of this is a high replication rate of cancer cells 
that needs to produce TGF-j) to sustain such growth rate. 

In our opinion the downregulation of the TGF-j) from the 
internalized TGF-j) has a defensive effect because it helps to 
reduce the internalization of TGF-j) [54]. In the bone niche, bone 
and cartilage contain large amounts of TGF-f) and target cells for 



TGF-j) activity. The autocrine and paracrine stimulation by TGF- 
ji is important for osteoblasts differentiation. The TGF-j) 
promotes osteoblast diflferentiation and the TGF-j) from the 
osteocytes in the bone make them lining up to the bone. In case of 
fracture the osteocytes secrete RANKL which activates the 
osteoclasts. The demolition of the bone around the fracture frees 
lot of TGF-j! which further attracts the osteoclasts but in the long 
term induces apoptosis. 

The model proposed suggests an ecological perspective of the 
cancer. Cancer cell changes are associated with alterations in the 
mechanical properties of the microenvironment; as tumor expands 
there is an increase in tissue compression and interstitial pressure, 
generating cell and tissue tension within the confined stroma. 
These forces induce the release and activation of various growth 
factors and subsequent changes of the contractility and viscoelas- 
ticity of tumor cells. In breast cancer, the different stages of cancer 
cell determine different mechanical interaction with basement 
membrane (BM) architecture and ECM. TGF-j) plays a key role 
in orchestrating the cell-ECM tension. A hyperplasia lesion 
typically involves the loss of normal cell polarization and 
organization, the changes in cell-cell contacts and cell-ECM 
interactions, which result in altered cellular tension and mechano- 
sensing and transduction. In carcinoma in situ lesions, cell polarity 
is lost and the lumen is filled by cells. This volume expansion and 
resistance from the BM and interstitial ECM lead to increased 
mechanical forces between tumor cells and the stromal matrix. 
Simultaneously, ECM components are abnormally deposited and 
remodeled, which results in increased ECM and tissue stiffness, 
and in turn, cell generated tension. In invasive lesions, tumor cells 
break down the BM and invade into the interstitial ECM. The 
reciprocal forces between tumor cells and the ECM continuously 
increase. The abnormal deposition and remodeUng of ECM 
collagen further increase ECM and tissue stiffness. Tumor cells 
generate greater tension in response to this increased mechanical 
stimulation. As tumor cells invade through the BM and ECM, they 
experience a range of different forces from the dense ECM 
network. 

Although TGF-j) is a growth inhibitor for most epithelial cells, it 
has multiple and often opposing effects depending on the tissue 
and the type of cells. Why TGF-j) is not regulated more tighfly? 
We believe that cancer is a disease related to the ageing; it is stiU 
very rare in young and mature organisms, while is very common in 
the elderly; so there is no great selection feedback and cancer is a 
mean to send in apoptosis an aged organism when most of the cells 
have accumulated mutations, so in some analogy it has the same 
role that TGF-j) has with single cells when energy becomes 
limited. 

Methods 

We have considered parameter estimates from experiments 
reported in literature and from published mathematical and 
computational models. Furthermore, some parameters are explic- 
itiy obtained from gene expression analysis, while some other are 
derived from the models so that their estimation as well as the 
models' outputs remain inside ranges of validity consistent with the 
biological phenomena. The values of the parameters and the 
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boundary conditions for the ODEs and DDEs used in the 
numerical simulations of the models are given in Tables (1, 2). 

In the cellular model, the values for the amount of the various 
isoforms of TGF-fS and receptors are extracted from gene 
expression data. In the specific case shown in Figure (3) and 
Table (2), where the cellular model and the tissue model are 
coupled together, we have considered the TGFB2 and TGFR3 
isoforms for the cytokine and its receptor respectively. 

The TGF-fi paracrine rate is derived as a direct conseriuence of 
the first Fick's law with an effective diffusion coefficient that takes 
into consideration the effects of the ECM. The order of magnitude 
for is estimated from the data in [55-58] relative to a healthy 
mammary epithelial cell. The reduction of volume in cancer 
mammary epithelial cells, as highlighted in [56], is the sign of cell 
cytoskeletal instabilities due to mutations. As consequence of the 
cell contractions, a small empty volume, which we have named 
diffusive region, appears between the cancer cell and the nearest 
cells of the tissue. The relation between the diffusive region volume 
and the neighbour volume is Fneigh « 7( I^cell + ^df) and Fceii is the 
averaged volume of a healthy mammary epithelial cell [56]. For 
the numerical values of ^syn linking the production of TGF-p to 
their gene expression, of the TGF-ji signaling rate to initiate the 
reactions with the Smad proteins in the cytoplasm and of the 
TGF-fi binding rate occurring at the cell surface, we refer to the 
experimented results in [16] with the assumption that Smad 
concentrations remain constant. The order of magnitude of 
mammary epithehal cell proliferation rate is based on the data 
reported in [59,60]. For the tissue model, the mutation rate and 
the probability, pt^, that a mammary cell undergoes malignant 
mutation showing a strong resistance to TGF-fS apoptotic 
signaling, are extracted from [21]. Based on the works [20,21], 
we have estimated the proliferation and the apoptosis exponents 
such that the local cell sub-population with highly malignant 
phenotype (<^ = 3) reaches the 10% of the local collective cell 
population approximatively after 2 years and the cell with 
phenotype (^ = 2 represents those cells with sufficient driver 
mutations so that the signaling induced by the TGF-P begin to 
change from anti-tumorigenic to pro-tumorigenic. Most of the 
parameters used in the bone niche model are taken from the 
following works: [47,53,61]. For the time immature osteoblasts 
take to differentiate into mature osteoblasts, the rate of osteocytes 
formation, the apoptosis rate of the osteoblasts and the rate at 
which the RANKL is released by the osteoblasts, we refer to the 
results in [46,62]. For the quantity of TGF-P stored in the bone, 
we refer to the work ofjanssens [63]. 

Gene Expression Data 

Despite the wealth of molecular data (such as sequence and 
gene expression data) and physiological and pathological data 
from different populations, there is a lack of cell abundance 
estimate in the different tissues. Therefore we extracted informa- 
tion from gene expression data. Gene expression patterns supply 
insight into complex biological networks. Gene expression 
profiling of the tumor microenvironment during breast cancer 
progression distinguishes breast carcinomas from normal breast 
tissues. We have re-analysed several gene expression data related 
to breast cancer dynamics from the Gene Expression Omnibus 
(http:/ /www.ncbi. nlm.nih.gov/geo/); after an exploratory analysis 
and literature analysis, we focused on the following datasets: 
(accession numbers): GSE14548, GSE33450 and GSE8977. These 
datasets originate from experimental design on early stages breast 
cancer progression and tumor microenvironment. Normalization 
procedures and statistical analysis are performed by using 
Bioconductor R packages [64]; the background correction and 



normalization is performed by using PLIER algorithm. PLIER 
algorithm produces an improved gene expression value [65] as 
compared to the other algorithms. It accomplishes this by 
incorporating experimental observations of feature behaviour 
Specifically, it uses a probe affinity parameter, which represents 
the strength of a signal produced at a specific concentration for a 
given probe. The probe affinities are calculated using data across 
arrays. The Bioconductor package limma was also used to 
calculate average expression levels, log fold changes and adjusted 
p-values for each probe. Standard anova and Box plots 
representation were used to analyse and check out visually the 
expression levels of these genes for different conditions. We have 
used gene expression averaged quantities to better unveil the 
functions of the TGF-ji in the cancer dynamics; nevertheless, the 
model proposed should be used with a pinch of salt. Single patient 
gene expression values may not be sufficient to catch the specificity 
of the evolution of the disease for each patient and the role of the 
TGF-P. Comorbidities, patients' previous medical conditions and 
the exact initial time of breast cancer formation may be unknown, 
but they affect the prediction of cancer evolution. We have also 
considered Transforming Growth Factor beta 3 (TGFB3) involved 
in cell differentiation and which interacts with TGF-fi receptor 2 
(TGFBR2), a tumor suppressor gene. 

Sensitivity Analysis 

The correctness of the proposed multi scale model and its 
predictive abilit)' on development and progress of the disease 
strongly depend on the knowledge of the parameters and the 
respective errors; hence, it is important to righdy identify and 
estimate the parameters whose the multi scale model is more 
sensitive. Using the Local Sensitivity Analysis method (LSA) in 
[66] , we compute the first order sensitivity indexes for each of the 
three models separately. Even though the non-linearity of the 
models suggests that the second order sensitivity indexes analysis 
could be relc\ ant to identify the most sensitive model's parameters 
couples, for sake of simplicity, we neglect the second order 
variance decomposition, and we show the LSA for all single 
parameters. The first order sensitivity index S, for a parameter Xi 
given the parameters X and the output Y = /(X) is a pure number 
defined as 5, = V{E(Y\Xi))/ V{Y), where E(Y\Xi) is the condi- 
tional expectation value obtained by sampKng all the parameters 
X from their prior distributions except Xi which is maintained 
constant, V(E(Y\Xi)) is the variance of E{Y\Xi) due to the 
variation of and V{Y) is the variance due to variation of aJl the 
parameters X. 

In the cellular model, the concentration of active TGF-P 
present in the extracellular region is strongly sensitive only to those 
parameters directly related to their gene expressions, Figure (5 .A, 
5.B). The same holds true for the TGF-f} receptors on the cell 
membrane. Figure (5.C, 5.D). Consequentiy, accurate information 
regarding the gene expression are much more important than 
other parameters to finely tune the outputs of the model and 
reduce the error on the capability of the model to be predictive. 
The fact that gene expressions play a fundamental role on the 
model is because they appear as terms of source of the proteins. 
Indeed, the other parameter that partially influence the model is 
the rate of synthesis, which also appears as factor in the source 
terms of equations (2, 3). Nevertheless, being the prior of the 
parameters uniformly distributed with a variance of 10% of their 
values expressed in Table (2), the sensitivity of the system to 

^syn tS 

approximatively two orders of magnitude less than to the gene 
expressions. On the other hand, the internalized compound not 
only is sensible to the gene expressions of the TGF-fi and its 
receptors, but also to the possibility cells have to internalize the 
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Figure 5. Cellular model LSA. The first order sensitivity index of the cellular model's parameters for (A,B) the TGF-fi (C,D) the receptors, (E,F) the 
TGF-fS internalized and (G) active TGF-fi in the diffusive region. On the left column the phenotype i)fi = 0 and on the right column ^ = 4 with 
exception of the TGF-ji in the diffusive region (G) vj/hich is independent of the phenotype. 
doi:10.1371/journal.pone.0088533.g005 
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Figure 6. Tissue model LSA. The first order sensitivity index of the tissue model's parameters for the tissue sub-population: (A) 1^ = 0, (B) ^ = 1, (C) 

1^ = 2, and (D) <^ = 3. 
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active form of TGF-fi through the receptors. In Figure (5.E, 5.F), 
we see tliat the quantity of TGF-ji accumulated in the cell 
cytoplasm depends on the capability of each specific cellular sub- 
population to respond to the TGF-ji signaling as well as on the 
capability of healthy cells to ubiquitinate the TGF-ji. The TGF-ji 
in the diffusive region depends on how fast it diffuses in the 
surrounding space and on its production by the healthy cells, 
Figure (5.G). Indeed, we are observing the system, which is in an 
out-of-equHibrium condition, at a time when the quantity of 
healthy cells is the predominant cell sub-population component of 
the total cell population; therefore, in this case, there are no 
contradictions in saying that the major sources of TGF-ji are the 
healthy cells. 

The sensitivity analysis for the tissue model presents a more 
complex scenario of the relevant parameters affecting the densities 
of each cell sub-population based on the cells phenotypes, see 
Figure (6). Due to the interaction driven by the paracrine 
exchange of TGF-P, the volume exclusion constraints between 
the cellular capacities and the strong non-linearity produced by the 
different cell capabilities of sensing/ responding to the TGF-fS, we 
see that the dynamical evolution of a cell sub-population with 
specific phenotype depends on parameters strictly non-related to 
that specific sub-population phenotype. The simplest case is given 
by the healthy cell sub-population which is prevalently sensitive to 
its proliferation/apoptosis rates, its exponential indexes, and the 
healthy cell capacity. Figure (6.A). Depending on the time at which 
the system is observed, we can see a more or less dependency of 
the healthy system on the proliferation rate of the pre-tumoral 
cells. The cell sub-populations characterized by driver mutations 
show all the same behaviours of being sensitive to the prolifera- 



tion/apoptosis rates and exponential indexes not only of their 
same phenotype, but also to those with less aggressive phenotypes. 
Figure (6.B, 6.C, 6.D). The small dependency on the parameters 
with phenotype = 2 is due to the particular transition role of this 
cell sub-population in respect to its response/sensing of the TGF-ji 
and consequently to the small range of the proliferation/ apoptosis 
indexes (meaning small variance of the prior distribution adopted 
in the LSA). All cell sub-populations do not show any sensitivity on 
the capacity C3, exception done for the highly aggressive tumoral 
cells which depend on both values of the cell capacities. 
Figure (6.D). 

From the LSA, Figure (6), we see the emersion of an explicit 
pattern on the parameters which reflects the behavioural evolution 
of the tissue in terms of dynamics. In fact, even though all the cells 
compete for limited room by obstructing or killing the cells with 
difiFerent phenotypes, the number of tumoral cells and their 
survival, during the early stage of cancer, are strongly linked to 
their phenotypical ancestors. Hence, the two major sub-systems of 
the tissue, healthy tissue and tumor, show the dual behaviour of 
two resources competing systems like in the predator-prey model, 
on one side and of two shearing TGF-ji systems in which one try 
to stabihze via TGF-ji the other sub-system without any success 
because the driver mutations in the aggressive cancer cells have 
changed the apoptotic signaling of the TGF-ji into a tumor 
promoter, on the other side. 

The LSA in the bone model shows two groups of variables: 
those which are mostly sensitive to the number of osteocytes in tiie 
BMU before a fracture occurs, but do not show any relevant 
sensitivity to the TGF-ji related parameters, and those variables 
which, even though depend on the osteocytes, are sensitive to the 
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Figure 7. Bone model LSA. The first order sensitivity indexes of the bone model's parameters for (A) osteocytes, (B) osteoclasts, (C) lining cells, (D) 
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delay induced by the TGF-fi as well as to the other parameters 
responsible to the high concentration of TGF-fl in the bone niche, 
Figure (7). The number of lining cells depends on the time they 
take to differentiate and hence they are sensitive to the TGF-fi, 
Figure (7.C). On the other hand, mature osteoblasts are not very 
sensitive to the TGF-j}, but are sensitive to the apoptotic rate of 
the lining cells because soon or later all the lining cell, which do 
not undergo apoptosis, foUow through maturation. Figure (7.D). 
We have also observed that the RANKL presents a strong 
sensitivity to the TGF-p which is an indication that the small 
quantity of RANKL released by the osteoblast compared to that 
released by the damaged osteocytes network is important to affect 
the dynamics of the bone remodelling, Figure (7.E). Consequently, 
the dependence of osteoclasts on TGF-ji is due to the RANKL, 
Figure (7.B). The mineral bone density is the variable with the 



highest sensitivity to the TGF-ji, Figure (7.F), showing that the 
resorption and remineralization are carefully synchronized by the 
TGF-P which regulates the time interval between osteoclasts and 
osteoblasts actions. 

It is important to highlight that the roles of osteoclasts and 
breast cancer cells, in the bone niche, are the same in terms of 
TGF-fi, but are extremely different in their resultant actions. The 
first (osteoclasts) produce TGF-fi to induce a delay necessary to 
complete the bone resorption and to carefully balance the 
maturation of lining cells, but at the same time, they are miners 
releasing a strong chemoattractant for the tumoral cells from the 
breast lobular duct. The last (breast cancer cells) remaining in the 
bone niche because rich of TGF-fi, which is a fundamental 
resource for their survival, release TGF-fi to prolongates the 
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mining of tlie osteoclasts, to make themselves space and to 
maintain their cell cycle progression. 

As pr(-\ i()u,sly stated the LSA is a useful tool to determine the 
importance of the parameters based on the fact that a parameter 
of the model with high sensitivity and low variance should be 
carefully chosen and/or measured because highly affecting the 
outcomes of the model while a parameter with low sensitivity and 
high variance does not influence the system dynamics. Further- 
more, we have used the LSA to reveal the relationships of and the 
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